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Abstract: We demonstrate an enhancement of the plane wave expansion 
method treating two-dimensional photonic crystals by applying Fourier 
factorization with generally elliptic polarization bases. By studying three 
examples of periodically arranged cylindrical elements, we compare our 
Q I approach to the classical Ho method in which the permittivity function is 

simply expanded without changing coordinates, and to the normal vector 
Qh| method using a normal-tangential polarization transform. The compared 

O . calculations clearly show that our approach yields the best convergence 

C^ \ properties owing to the complete continuity of our distribution of polariza- 

tion bases. The presented methodology enables us to study more general 



^ I systems such as periodic elements with an arbitrary cross-section or devices 

^ ^. such as photonic crystal waveguides. 
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1. Introduction 

Photonic crystals (PhCs) are modern artificially structured systems of a broad interest from 
many viewpoints of fundamental research and applications d El [3] IH . Studying their electro- 
magnetic properties is significantly important for developing PhC-based promising applications 
such as purely optical integrated circuits |[ll[6l[7l, artificial metamaterials with high tunabil- 
ity |[8l|9j[T0l[TTl, high-sensitivity PhC biosensors |[l2l[T3l, or devices based on phenomena not 
accessible in conventional media ||T4l [15] [161 • It is also necessary for accounting for structural 
colors of wings of butterflies or beetles, feathers of birds, or iridescent plants |[T7l[T8l . 

Owing to the PhC periodicity, the plane- wave expansion (PWE) method is conventional for 
calculating PhC modes and photonic band structures. Its essence is the Fourier expansion of 
electromagnetic fields and material parameters, which is not only a counterpart of the PWE 
method for electronic crystals, but it also advantageously follows the classical coupled-wave 
theory developed for diffraction gratings during last several decades |[T9| |20l |2T1 . One of the 
most important discoveries of the grating theory are the Fourier factorization (FF) rules for 
expanding the ratios of the discontinuous permittivity and the discontinuous electric field [22J . 
The three FF theorems were first formulated for one-dimensional (ID) gratings and then gen- 
eralized to two-dimensional (2D) systems ||23]| , arbitrary periodic reliefs ll24l . anisotropic ||25]| 
and slanted ll26l periodic structures, their various combinations ||27l |28l |29l, and other ap- 
plications Oni |3T| |32l . As will be demonstrated in this paper, the FF rules establish suitable 
principles for an accurate solution of Maxwell's equations in the plane- wave basis, which con- 
siderably enhances the numerical performance compared to previously used implementations. 

In the case of 2D-periodic structures, the FF rules were first applied to "zigzag" Fourier 
expansions ll23ll , which yielded an improvement for rectangular dots or holes, but were not 
sufficient for the staircase approximation of round elements. After that a coordinate transform 
was applied to treat individually the normal and tangential components of the electric field on 
ID sinusoidal-relief gratings, which enabled the application of the correct rule for each field 
component ll24l . Later David et al. ll33]| utilized the normal-tangential field separation to 2D 
PhCs composed of circular and elliptical holes. Similarly, Schuster and colleagues ll34l ap- 
plied this method to 2D gratings, and also suggested more general distributions of polarization 
bases ll35]| . These approaches, always dealing with linear polarizations, enabled a significant 
improvement of the convergence properties, but ignored the fact that the transformation matrix 
between the Cartesian and the normal-tangential component bases of polarization became dis- 
continuous at the center and along the boundaries of the periodic cell, which slowed down the 
resulting convergence. To overcome these discontinuities, a distribution of more complex (i.e., 
generally elliptic) polarization bases was recently suggested to improve optical simulations of 
2D gratings |36|. 

Here we are dealing with 2D PhCs made as 2D-periodic elements of the circular cross- 
section, for which we implement the PWE method by using FF with complex polarization 
bases. From the mathematical point of view this complexity only requires complex- valued 
Jones vectors, so that the generalization from the previous method of David is surprisingly 
simple. Both methods are compared with respect to their convergence properties, together with 
the classical method of Ho [3 1 (where the electric permittivity is expanded without any coordi- 
nate transform), to show that the method presented here yields the best performance. 

2. PWE method for 2D structures 

We study a 2D PhC composed of infinite cylinders with a circular cross-section with either 
square [Fig. [Ha)] or hexagonal [Fig.[nb)] periodicity. While the structure is invariant along the 
z-axis, the periodicity is assumed in the x- and j-axes directions. For the square symmetry the 
unit cell has the dimensions ax = ay = a, and for the hexagonal symmetry it can be chosen as an 
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Fig. 1. Two studied configurations of 2D PhCs, with square (a) and hexagonal (b) sym- 
metry, with denoted square and rectangular unit cells, respectively. The corresponding first 
Brillouin zones of the reciprocal space are depicted in (c) and (d), respectively, where 
the emphasized triangles denote irreducible fractions. Note that for calculations we use a 
rectangular unit supercell (solid rectangle) rather than a hexagonal primitive cell (dashed 
hexagon) in the case of the hexagonal lattice (b) so that the corresponding first Brillouin 
zone is the solid rectangle instead of the dashed hexagon in (d), which creates a folded band 
structure, from where we choose values corresponding to the usual convention. 

a-by-a>/3 rectangle. The corresponding first Brillouin zones of the reciprocal space are depicted 
in Figs.IHc) andlljd), together with the symmetry points F, X, M, and F, M, K, respectively. 
According to this 2D-periodic structure we define and expand the relative permittivity function 
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where e^„ are its Fourier coefficients and G^ = 27t/ax and Gy = iTc/uy are the reciprocal lattice 
vectors. 

Maxwell's equations for the //-polarization (where the electric field E is in the x-y plane 
whereas the magnetic field H is along the z-axis) are 
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where CO is the photon frequency, and £o and jiq are the permittivity and permeability of vacuum, 
or 
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Note that here we do not study the ^^-polarization (for which the electric field E is along the 
z-axis) because in that case the Ho method satisfies the correct Fourier factorization rules. 

Assuming an in-plane anisotropy of the relative permittivity function, defining a scaled elec- 
trical displacement D, 
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defining a 2-by-2 matrix of electrical impermittivity TJ 
Eq. ([5]) yields 



e ^ and substituting Eq. Q into 
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Fig. 2. Schematic description of the polarization distributions in Models A (a), B (b), C (c), 
and C (d) within the first quadrant of the periodic cell, as functions of the polar coordinates 
re^^ = x-\- iy. In (a) the x-y Cartesian basis is uniform; in (b) the polarization vectors are 
normal (u) and tangential (v) to the cylindrical element and constant along lines coming 
from the center; in (c) and (d) the u, v polarizations are in general elliptic which enables 
their continuity. 



where r\jk are the components of the electrical impermittivity, c = (£oMo)~^^^ is the light speed 
in vacuum, and the factor (O^ /c^ is calculated as the eigenvalues of the operator on the left-hand 
side of Eq. ©. 

According to the Floquet-Bloch theorem, in the reciprocal space the magnetic field under- 
goes the Fourier expansion 
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where pm = kx-\-mGx, qn = ky-\-nGy. Analogously, the Fourier components of E^, Ey, D^, by 
are denoted ^x,mw, ey^mn, dx,mn, dy^mn- Then, as demonstrated in Appendix, the operators of the 
partial derivatives dx, dy become the diagonal matrices — /p, — /q, the operators of multiplication 
by functions such as £xx or rfxx become matrices denoted [[e^x]] or [[t7xx]] , and the field compo- 
nents Ex, Ey, Hz become column vectors denoted [Ex], [Ey], [H^. Using these definitions, the 
multiplications in Eq. (|7]) can be treated either by the Laurent FF rule, [e^x^x] = [[^xx]] [^x] etc., 
or by the inverse FF rule, [e^x^'x] = [[V^xx]]"^ [Ex]- 

3. Methods of Fourier factorization 

3.1. Elementary (Cartesian) method (Model A) 

In this article we compare several models corresponding to different FF approaches. First, 
Model A assumes the solution in the basis of the x and y polarizations uniform within the 
periodic cell [Fig.[2ta)], where in accordance with the Ho method |3] we choose the Laurent 
rules 

[Dx] = [eEx] = [[e]][Exl (10) 

[by] = [eEy] = [[E]][Ey]. (11) 

The components of the electric impermittivity in Eq. ([5]) then becomes [[e]]~^ for the cases of 
^xx, ^yy^ and zero for the cases of r{xy^ Vyx^ or 
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In Fig. 12 a) we depict the uniform polarization basis within the first quadrant of the square 
periodic cell, < x < £>(0), < j < D{n/2), where £>((/») is the distance between the cell's 



center and its boundary (a function of the polar coordinate (j)). We also denote R the radius of 
the cylindrical element. 

3. 2. Normal vector method (Model B) 

According to Li ll22l the Laurent rule is valid (for the reason of uniform convergence) for 
multiplying functions that possess no concurrent discontinuities, whereas the inverse rule is 
used for functions whose product is continuous. Obviously, neither the Laurent rule nor the 
inverse rule is correct for both products in Eqs. (fTOl) and (fTTt . because both pairs of functions 
have concurrent discontinuities and both products Dj and Dy are discontinuous as well. On 
the other hand, by an appropriate change of the polarization bases at all points (using a space- 
dependent Jones matrix transform F), 
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we can treat independently the normal (u) and tangential (v) components of the fields by the 
correct rules. 
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The field components £„, Du are normal to the discontinuities of the relative permittivity func- 
tion, while Ey, Dy are tangential. The FF rules used in Eqs. (fT4l) and ([TSll are justified simply 
because Ey and D„ are continuous. 

A suitable distribution of the matrix F within the periodic cell can obviously be the rota- 
tion (SS 

cos (/» — sin (^ 
sin (j) cos (j) 



(16) 



where the polar angle 0(x,};) is in the first cell distributed according to the polar coordinates 
re^^ =x-\- iy, and then periodically repeated over the entire 2D space. This enables defining the 
matrices [[c]], [[i-]] from the corresponding 2D-periodic functions c = cos0, ^ = sin0. 

Let u and v be the two columns of the matrix F, both being mutually orthogonal basis vec- 
tors of linear polarization 137]. From the above definitions we see that u is a polarization vector 
normal to the structure discontinuities, whereas v is tangential. In Fig. [2tb) we depict the dis- 
tribution of the u-v basis within the first quadrant of the periodic cell; the polarization vectors 
are constant along the lines of the constant azimuth (0 = const) and rotate as (j) increases. The 
distribution of the rotation of u within the square periodic cell is displayed in Fig. Ha), from 
where it is obvious that the matrix function F(x, j) has no discontinuities concurrent with the 
electric field, so that we can use both Laurent and inverse rules for the transformation of polar- 
ization, e.g.. 
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Combining Eqs. ([T4b . (fTSl) , and (fTTl) yields 
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Fig. 3. Distribution of the rotation and ellipticity of the basis polarization vector u for 
the presented models. The rotation for Models B and C is in (a) [(e) for the hexagonal 
periodicity]; the ellipticity for Model C is in (b) [(f) for the hexagonal periodicity]; and the 
rotation and ellipticity for Model C are in (c) and (d), respectively. The color scales for both 
the rotation and the ellipticity are on the right (in degrees). The structure discontinuities (the 
circular boundaries of the periodic elements) are plotted as black or white circles. Notice 
that u is always linear along and normal to the circle. 



from where we derive the electric impermittivity in the reciprocal space (corresponding to 
IModel B) 
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whose components are immediately applicable to Eq. ([8]). 



3.3. Method with elliptical polarization bases (Model C) 

The above approach (IModel B) only deals with linear polarizations and thus suffers from the 
fact that the matrix function F(x, j) has a singularity at the center of the periodic cell and other 
discontinuities along the cell boundaries. This slows down the convergence of the numerical 
implementation, as will be evidenced below. 

As suggested in the previous work for the case of diffraction by 2D gratings 1361 , we can 
make F continuous by using complex functions ^ and ^ or, in other words, by defining u, v as 
complex vectors corresponding to generally elliptic polarizations. 






-c* 



(21) 



(which are still orthogonal). By means of rotation 6 and ellipticity (o we define the first basis 
vector 
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where 
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r f(l+cosf) {r<R) 



Here R denotes the radius of the circular element and 
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is the distance from the cell's center to its edge. In Eq. (l22l) the Jones vector on the right rep- 
resents a polarization ellipse (with ellipticity (o) oriented along the x coordinate, the matrix in 
the middle rotates this polarization by the azimuth 0, and the factor e^^ preserves the continuity 
of the phase at the center and along the boundaries of the cell. This continuity can be easily 
checked by evaluating the limits 
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which is the vector of left circular polarization (independent of (j)). 

The distribution of the u-v polarization basis within the first quadrant of the periodic cell 
is schematically depicted in Fig. [2c), where the azimuth of the polarization ellipse is constant 
along the lines coming from the cell's center, which is similar to Model B; hence. Fig. Oa) 
serves for both Model B and Model C. However, the ellipticity is now zero (corresponding 
to linear polarization) only on the boundaries of the circular element, has the maximum value 
(7l/4 for circular polarization) at the cell's center and along its boundaries, and continuously 
varies (with a smooth sine dependence) in the intermediate ranges [Fig. Ob)]. Finally we obtain 
a smooth and completely continuous matrix function F(x,j), which is analogously used to 
calculate the impermittivity in the reciprocal space 
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In the case of the hexagonal periodicity we define u and the other periodic quantities inside 
one hexagon (half the area of the rectangular unit cell) where we can use formally the same 
equations as above, except for 

max cos (O — -^)\ 

which is now the distance from the hexagon's center to its edge. Here a is the hexagon's shortest 
width (equal to the width Gx of the rectangular cell). After periodic repeating over the entire 2D 
space, the azimuth and ellipticity distributions of the basis vector u become those displayed 
in Figs. Oe) and Of), respectively. Note that although the color distribution in Fig. Oe) is 
discontinuous along the hexagon's boundaries, the vector u is actually continuous, because the 
color difference is only due to the 180° change. Moreover, the vector u in the limit r — > D{(l)) 
corresponds to left circular polarization, so that its azimuth becomes irrelevant. 

3.4. Modified method for densely arranged elements (Model C) 

To analyze a more complicated situation, we consider a PhC with square periodicity where 
circular elements are densely arranged near each other, i.e., where the radius R is almost the 
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Fig. 4. Amplitude distribution of the scaled magnetic field \Hz{x,y)\ within one cell of the 
real space for Samples SI (a), S2 (b) and H (c). The color scale, same for each sample, 
is in the top left corner of (a). The subfigures represent eigenmodes distinguished by the 
symmetry point letters (F, X, M) and band numbers (1-4). The structure discontinuities 
(boundaries of the rods or holes) are plotted as white or black circles. 



half width a/2 of the periodic cell. Then the convergence properties of F becomes worse, which 
affects all the derived quantities. For this reason we again redefine the polarization distribution, 
as suggested in Fig. [2d). For the modified Model C we define u to be still same inside the 
circle (r < R), but different outside. Assuming the rotation and ellipticity along the boundary 
of the square cell 
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(where "round" denotes rounding towards the nearest integer), we define the rotation and el- 
lipticity outside the circle (r > R) as 
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Assuming otherwise the same Eqs. ([T9l ), (|2T1) , (l22l) , and (|25]) , we obtain for [[tj]]c/ formally the 
same matrix as in Eq. (|27]) , except that the functions t, and ^ are now derived from different 
azimuth and ellipticity distributions of u [displayed in Figs.Oc) and^d), respectively]. Note 
that u is again continuous along the cell's boundaries; to evaluate its precise limits [when x -^ 
±Z)(0), y = const or j ^ ±D{(j)/2), x = const] is now more complicated. 

4. Numerical examples 

4.1. Description of simulations 

We examine the numerical performances of all the presented models on three samples of 2D 
PhCs, for which we calculate the eigenfrequencies co^ (where the band number K = 1 stands 
for the lowest eigenfrequency, K* = 2 for the second lowest, etc.) and the corresponding eigen- 
vectors [Hz] K of Eq. ([8]). All convergences will be presented according to the maximum Fourier 
harmonics retained inside the periodic medium (M, N defined in Appendix), which will be kept 
same for the x and y directions (M = N). 

First, Sample SI is a square array of cylindrical rods of the circular cross-section with the 
diameter 2R = 500 nm, square period a = 1000 nm, relative permittivity of the rods £i = 9, and 
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Fig. 5. Convergence properties of normalized eigenfrequencies (coa/lnc) calculated for 
selected bands of Sample SI. The modes r2, FS are on the left part of the figure; XI, X4 in 
the middle; and M2, M4 on the right. The Models A, B, and C are compared to each other, 
plotted as crosses, squares, and circles, respectively. 



relative permittivity of the surrounding medium corresponding to vacuum (£2 = 1). For clarity 
we display the first several eigenmodes at the points of symmetry F, X, M in Fig. HJa). The 
color scale of the modes corresponds to the scaled amplitude distribution of the magnetic field 
\Hz{x^y)\ within one cell of the real space (calculated by Model C with M = N = 12). Each 
eigenmode is denoted by the symmetry point letter and the band number (e.g., X2 denotes the 
second band at the point X). By comparing the calculated eigenfrequencies and the amplitude 
distributions we find that the pairs of eigenmodes F3-F4 and M2-M3 constitute frequency 
doublets, i.e., the frequencies for one pair are very close to each other (differences are due to 
numerical errors). Note that although the amplitudes \Hz{x,y)\ of the two eigenmodes within a 
doublet seem identical [e.g., the doublet F3-F4 in Fig.JHa)], the two actual eigenmodes, which 
are complex- valued distributions Hzix^y), are indeed different, being two linearly independent 
eigenvectors of a degenerate eigenvalue. 

Similarly, for Sample S2 we assume exactly the same parameters except the diameter of the 
rods, now being 2R = 900 nm. This corresponds to densely arranged elements (the distance 
between two adjacent rods is only 100 nm). The amplitude distributions for the same bands at 
the same symmetry points are displayed in Fig.|Hb) (calculated by Model C with M = A/^ = 12), 
where we again observe the frequency doublets in the same bands. 

Finally, for Sample H we consider a hexagonal array of cylindrical holes of the circular cross- 
section with the diameter 2R = 600 nm, hexagonal periodicity a = 1000 nm (corresponding to 
the rectangular cell of the dimensions aj = 1 fim, ay = VS /im), relative permittivity of the 
holes corresponding to vacuum (£i = 1), and relative permittivity of the substrate medium 
(surrounding holes) £2 = 12. Since the lowest (nonzero) eigenfrequencies at the F and K points 
are quite high, we confine ourselves only to the first three eigenmodes at the point M, whose 
amplitude distributions (within the rectangular cell) are displayed in Fig.|4tc). 




Fig. 6. Convergence properties of normalized eigenfrequencies (coa/lnc) calculated for 
the modes Ml (left), M2 (middle), and M3 (right) of Sample H. The Models A, B, and C 
are plotted as crosses, squares, and circles, respectively. 



4.2. Results and discussion 

We carefully compared the numerical effectivity of the presented models for all the defined 
samples, and below we present the convergence results for selected bands. 

For Sample SI the difference among the performances of Models A, B, and C was found as 
follows: The values of Model B always converge much faster than those of Model A, which is a 
result expected from previous reports 1331 . The values of Model C converge significantly faster 
than those of Model B for the modes r3, r4, XI, X4, and M2-M4. This can be explained by 
the discontinuities of the matrix function ^{x^y) along the boundaries and at the center of the 
periodic cell, which are responsible for the Gibbs phenomenon appearing in the partial Fourier 
sums in Eq. ([2Qb . On the other hand, both Models B and C exhibit similar convergence prop- 
erties for the modes r2, X2, X3, and Ml, where the amplitude distribution has almost circular 
symmetry and is nearly zero near the cell's boundary [Fig.|4l;a)], so that the discontinuities do 
not manifest themselves. 

A few examples of compared convergences are displayed in Fig. [5l where a considerable 
improvement owing to Model C is visible everywhere except for the T2 mode. Nevertheless, 
according to the scale of the frequency axis at the mode T2 no improvement by Model C is 
necessary since Model B already yields the 5-digits precision. For the modes r3, XI, and X4 
the precision is improved by at least one order. 

For Sample H the convergence of Model A is also much slower than those of Models B 
and C. In Fig. [6] the values of Model A are even beyond the displayed range of the graphs for 
the modes M2 and M3. Model C yields the best results for the modes Ml and M2. On the other 
hand. Models B and C converge similarly for the mode M3; obviously a higher range of N 
would be required to compare the convergences more adequately. 

All the eigenmodes r2-r4, X1-X4, and M1-M4 of Sample S2 were carefully analyzed with 
respect to the convergence properties of Models A, B, C, and C. Here the values of Model C, 
as expected, converge fastest for all the modes except Ml, whereas Model A is again the worst. 
The eigenmode of the mode Ml is again circularly symmetric and has negligible values near 
the cell's boundaries [Fig.jljb)], so that Model B is not affected by the discontinuities of the 
polarization transform and yields precise values. On the other hand. Model C is found inap- 
propriate for this sample (yielding even worse convergence than Model B) because the rapid 
variation of the ellipticity of u near the cell boundaries between two adjacent elements (which 
are very close to each other) requires more Fourier components than the weak discontinuity of 
the linear polarization u in Model B. This problem is removed in Model C, as obvious from 
Figs.[2l;d)and[ad). 
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Fig. 7. Convergence properties of normalized eigenfrequencies {(Oa/lnc) calculated for the 
modes XI (left), X2 (middle), and X3 (right) of Sample S2. The Models A, B, C, and C 
are plotted as crosses, squares, circles, and triangles, respectively. 



Three examples of convergences (for the modes X1-X3) are shown in Fig. [71 which clearly 
justifies the suitability of Model C. Notice again that the values of Model A are in two cases 
(modes XI and X2) beyond the displayed range of the graphs (with error higher by about two 
orders). In the case of the mode X3 the values of Model A do not differ so much, which can 
be explained by nearly zero values of the eigenmode X3 at all points of PhC discontinuities, 
which is obviously not the case of the modes XI and X2 [Fig.jljb)]. 



5. Conclusions 

We have derived three models with respect to different methods of applying the FF rules and 
compared their numerical performances. Model A (equal to the Ho method) exhibits the worst 
convergence properties. Model B gives the best performance in a few cases where the eigen- 
mode is circularly symmetric and has a negligible amplitude near the cell's boundaries (then the 
problems coupled with polarization discontinuities vanish in zero fields). Model C was found 
to converge significantly faster than Model B for PhCs without dense arrangement of elements. 
On the other hand, for dense PhCs, where Model C was found inappropriate, the best results 
were yielded by Model C, which combines the best properties of Models B and C. We can 
conclude that the method of FF with complex polarization bases, represented here by Models C 
and C, produces an enhancement of the PWE method, provided that we choose an appropriate 
distribution of the polarization basis according to the sample under study. 

The method can be straightforwardly generalized to PhCs composed of elements of an arbi- 
trary cross-section. It is possible by replacing the constant value of the element's radius R with 
a function 7?(0), and by a further generalization of the azimuth 0(r, 0) and ellipticity ^(r, 0) of 
the vector u to make it again normal to and linear at the element boundaries (and continuous 
in the entire 2D space). Moreover, the method can also be used to PhC-based devices such as 
PhC waveguides by applying the demonstrated methodology to the device supercell (similarly 
as we were here dealing with the rectangular cell of the hexagonal PhC). It is particularly ad- 
vantageous for devices where high accuracy is required, e.g., for analyzing defect modes near 
photonic band edges |[38l[39l|40l|, and for large devices for which the available computer mem- 
ory enables calculations with only a few Fourier harmonics, e.g., PhC-based fibers with large 
cladding. 
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Appendix: Matrices for 2D structures 

Here we follow the matrix formation procedures and notation as explained by Visnovsky and 
Yasumoto BB . Let f{x^y) be one component of the pseudo-periodic electric field inside the 
medium described by the 2D-periodic function e(x, j). Their Fourier expansions are written 

e{x,y) = £ £m„e-'('"°-^+««^^), (33) 

m,n=—oo 

f{x,y) = £ fmne-'^P-'+'-yK (34) 



h{x,y) = 


V /? p-KPmX^qny) 




m,n=-oo 


gx{x,y) = 


+ 00 

V P ^-iiPmX^qny) 
L^ 8x,mn ^ 




m,n=—oo 


gyi^^y) = 


+ 00 

V P ^-KPmX^qny) 

L^ 6y,mn ^ 



We will here derive the matrix expressions of the fundamental relations 

h{x,y) = £{x,y)f{x,y), (35) 

gx{x,y) = dj(x,y), (36) 

gy{x,y) = dyf{x,y), (37) 

i.e., the relations of multiplication by a function and applying partial derivatives. Assuming the 
expansions of the new functions 

(38) 
(39) 
(40) 

we rewrite Eqs. ([35l)-([37b using the convolution rule, and applying the partial derivatives as 
follows: 

hmn = 2^ ^m-k,n-lfkU (41) 

k,l=-oo 

gx,mn ^ iPmJmn^ V^^) 

gy^mn ^ l^nlmn- \^^) 

Assuming furthermore a finite number of the retained Fourier coefficients, i.e., using the sum- 
mation Y.m=-M^n=-N' ^^ ^^^ rcnumbcr all the indices so that instead of a couple of two sets 
m G {— M, —M-\- 1, . . . , M} and n G {—N^ —N-\- 1, . . . , A/^} we have a single set of indices 
a G {1, 2, ... , amax}, with Omax = (2M+ 1)(2A^+ 1), related 

a{m,n) = m + M+ 1 + (^ + A^)(2M+ 1), (44) 

n{a) = (a-l)div(2M+l)-A^, (45) 

m(a) = (a-l)mod(2M+l)-M, (46) 



where "div" denotes the operation of integer division and "mod" the remainder (the modulo 
operation). Then we can rewrite Eqs. (l4T1) -(l43l) into the matrix relations 

[h] = Mifi (47) 

[gx] = -W], (48) 

[gy] = -mlf], (49) 

where [/], [h], [gx], and [gy] are column vectors whose ath elements are the Fourier [m^n] 
elements of the functions /, h, gj, and gy, indexed by a{m^n) defined in Eq. (l44l) , and where 
[[e]] , p, and q are matrices whose elements are defined 

Ma/3 = ^m{a)-m{(5),n{a)-n{(5)^ (50) 

Pap = Pm{a)Saf5i (51) 

Qa(5 = ^n{a)Sa(5^ (52) 

where the indices on the right hand parts are defined by Eqs. (l45l) and (l46l) and where So^^g 
denotes the Kronecker delta. As a summary we can say that the multiplication by a function 
is in the reciprocal space represented by the matrix [[e]] (in the sense of the limit Omax -^ °°), 
which is a generalization of the Toeplitz matrix used for ID periodicity, and that the partial 
derivatives are represented by the diagonal matrices — /p and — /q. 



